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Abstract. The aim of this work is to provide fast and accurate approx- 
imation schemes for the Monte Carlo pricing of derivatives in LIBOR 
market models. Standard methods can be applied to solve the stochas- 
tic differential equations of the successive LIBOR rates but the methods 
are generally slow. Our contribution is twofold. Firstly, we propose an 
alternative approximation scheme based on Picard iterations. This ap- 
proach is similar in accuracy to the Euler discretization, but with the 
feature that each rate is evolved independently of the other rates in the 
term structure. This enables simultaneous calculation of derivative prices 
pH , of different maturities using parallel computing. Secondly, the product 

Q) ■ terms occurring in the drift of a LIBOR market model driven by a jump 

process grow exponentially as a function of the number of rates, quickly 
rendering the model intractable. We reduce this growth from exponen- 
tial to quadratic using truncated expansions of the product terms. We 
Q^| include numerical illustrations of the accuracy and speed of our method 

pricing caplets, swaptions and forward rate agreements. 
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1. Introduction 
The LIBOR market model (LMM) has become a standard model for the 



> 

\o 

^-v . pricing of interest rate derivatives in recent years. The main advantage of 

. | this model in comparison to other approaches is that the evolution of dis- 

^^ ■ cretely compounded, market-observable forward rates is modeled directly 

^^ , and not deduced from the evolution of unobservable factors. Moreover, the 

log-normal LIBOR model is consistent with the market practice of pricing 
caps according to Black's formula (cf. Black 1976). However, despite its ap- 
parent popularity, the LIBOR market model has certain well-known pitfalls. 
5_i ■ On the one hand, the log-normal LIBOR model is driven by a Brownian 



motion, hence it cannot be calibrated adequately to the observed market 
data. An interest rate model is typically calibrated to the implied volatility 
surface from the cap market and the correlation structure of at-the-money 
swaptions. Several extensions of the LIBOR model have been proposed in 
the literature using jump-diffusions, Levy processes or general semimartin- 
gales as the driving motion (cf. e.g. Glasserman and Kou 2003, Eberlein and 
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Ozkan 2005, Jamshidian 1999), or incorporating stochastic volatility effects 
(cf. e.g. Andersen and Brotherton-Ratcliffe 2005). 

On the other hand, the dynamics of LIBOR rates are not tractable under 
forward measures due to the random terms that enter the dynamics of rates 
during the construction of the model. In particular, when the driving process 
has continuous paths the dynamics are tractable under their corresponding 
forward measure, but not under any other forward measure. When the driv- 
ing process is a general semimartingale, then the dynamics of LIBOR rates 
are not even tractable under their very own forward measure. Consequently: 

(1) if the driving process is a continuous semimartingale caplets can be 
priced in "closed form", but not swaptions or other multi-LIBOR 
derivatives; 

(2) if the driving process is a general semimartingale, then even caplets 
cannot be priced in closed form. 

The standard remedy to this problem is the so-called "frozen drift" ap- 
proximation; it was first proposed by Brace et al. 1997 for the pricing of 
swaptions and has been used by several authors ever since. Brace et al. 
(2001), Dun et al. (2001) and Schlbgl (2002) argue that freezing the drift 
is justified, since the deviation from the original equation is small in several 
measures. 

Although the frozen drift approximation is the simplest and most pop- 
ular solution, it is well-known that it does not yield acceptable results, es- 
pecially for exotic derivatives and longer horizons. Therefore, several other 
approximations have been developed in the literature. In one line of research 
Daniluk and Gatarek (2005) and Kurbanmuradov et al. (2002) are looking 
for the best lognormal approximation of the forward LIBOR dynamics; cf. 
also Schoenmakers (2005). Other authors have been using linear interpola- 
tions and predictor-corrector Monte Carlo methods to get a more accurate 
approximation of the drift term (cf. e.g. Hunter et al. 2001 and Glasserman 
and Zhao 2000). We refer the reader to Joshi and Stacey (2008) for a de- 
tailed overview of that literature, some new approximation schemes and 
numerical experiments. 

Although most of this literature focuses on the lognormal LIBOR market 
model, Glasserman and Merener (2003b, 2003a) have developed approxima- 
tion schemes for the pricing of caps and swaptions in jump-diffusion LIBOR 
market models. 

In this article we develop a general method for the approximation of the 
random terms that enter the drift of LIBOR models that is suitable for 
parallel computing. In particular, using Picard iterations we develop generic 
approximation schemes that decouple the dependence between LIBOR rates. 
Therefore individual rates in the tenor can be evolved independently in a 
Monte Carlo simulation. In addition, we treat a problem specific to LIBOR 
models with jumps; namely that the complexity of the drift term grows 
exponentially in the number of tenor dates. We expand and truncate the 
drift term, which yields a highly accurate approximation, while the gain in 
computational speed is immense. We illustrate the accuracy and efficiency 
of our method in an example where LIBOR rates are driven by a normal 
inverse Gaussian process. 
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The method we develop is universal and can be applied to any LIBOR 
model driven by a general semimartingale. However, we focus on the Levy 
LIBOR model as a characteristic example of a LIBOR model driven by a 
general semimartingale. 

The article is structured as follows: in section 2 we review time-inhomoge- 
neous Levy process, and in section 3 we revisit the Levy LIBOR model 
and explain in detail the computational problems. In section 4 we derive 
the Picard approximation scheme and the drift expansions. In section 5 
we briefly describe the main derivatives on LIBOR rates. Finally, section 6 
contains a numerical illustration. 

2. Levy processes 

Let (f2, F, F, P) be a complete stochastic basis, where F = Ft* and the 
filtration F = (Ft)t£[o,T„] satisfies the usual conditions; we assume that T* G 
Mj>o is a finite time horizon. The driving process H = (H t )o<t<T t is a process 
with independent increments and absolutely continuous characteristics; this 
is also called a time-inhomogeneous Levy process. That is, H is an adapted, 
cadlag, real- valued stochastic process with independent increments, starting 
from zero, where the law of Ht, t G [0, T*], is described by the characteristic 
function 

TE[e iuHt ] = exp I I [i6 s u - ju 2 + I ' (e iux - 1 - iux)F s (dx)\ ds J ; (2.1) 

\0 R / 

here bt G K, c± G M^o and Ft is a Levy measure, i.e. satisfies Ft({0}) = 
and J" K (1 A \x\ 2 )F t (dx) < oo, for all t G [OjT*]. In addition, the process H 
satisfies Assumptions (AC) and (EM) given below. 

Assumption (AC). The triplets (bt,ct,Ft) satisfy 

I (\h\ + (k+ f(lA\x\ 2 )F t {dx)\dt<oo. (2.2) 

K 

Assumption (EM). There exist constants M, e > such that for every 
u G [-(1 + e)M, (1 + e)M] =: M 

f I e ux F t (dx)dt<oo. (2.3) 

{|*|>1} 

Moreover, without loss of generality, we assume that fn x \ > i} e ux F t (dx) < oo 
for all t G [0, T*] and u G M. 

These assumptions render the process H = (H t )o<t<T t a special semi- 
martingale, therefore it has the canonical decomposition (cf. Jacod and 
Shiryaev 2003, II. 2. 38, and Eberlein et al. 2005) 

H = f b s ds + I ^c~dW s +!( x(fi H - u)(ds, dx), (2.4) 

OR 
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where fj, is the random measure of jumps of the process H, v is the P- 
compensator of fj, H , and W = (Wt)o<t<T* is a P-standard Brownian motion. 
The triplet of predictable characteristics of H with respect to the measure 
P, T(H\W) = (B,C,u), is 

B = b s ds, C = c s ds, i/([0,-] x A) = I I F s (dx)ds, (2.5) 

A 

where A G B(M); the triplet (b, c, F) represents the local characteristics of H. 
In addition, the triplet of predictable characteristics (B, C, v) determines the 
distribution of H, as the Levy-Khintchine formula (2.1) obviously dictates. 
We denote by k s the cumulant generating function associated to the in- 
finitely divisible distribution with Levy triplet (b s , c s ,F s ), i.e. for z G M and 

sG [0,7;] 

K,(z) := b s z + ^-z 2 + f(e zx - 1 - zx)F s (dx). (2.6) 

R 

Using Assumption (EM) we can extend k s to the complex domain C, for 
z G C with Kz G M, and the characteristic function of i7f can be written as 

t 
E [e-.]^xp(/„,(,„) d ,). (2.7) 



If i7 is a Levy process, i.e. time- homogeneous, then (b s , c s ,F s ) - and thus also 
k s - do not depend on s. In that case, k equals the cumulant (log-moment) 
generating function of H\. 

3. The Levy LIBOR model 

3.1. Model description. The Levy LIBOR model was developed by Eber- 
lein and Ozkan (2005) following the seminal articles on LIBOR market mod- 
els driven by Brownian motion by Sandmann et al. (1995), Miltersen et al. 
(1997) and Brace et al. (1997); see also Glasserman and Kou (2003) and 
Jamshidian (1999) for LIBOR market models driven by jump processes and 
general semimartingales respectively. The Levy LIBOR model is a market 
model where the forward LIBOR rate is modeled directly and is driven by 
a time-inhomogeneous Levy process. 

Let = To < T\ < ■ ■ ■ < T/v < Tjv+i = T* denote a discrete tenor 
structure where 5i = Tj + i — Tj, i G {0, 1, . . . , N}. Consider a complete sto- 
chastic basis (f2, F, F,P^) and a time-inhomogeneous Levy process H = 
(Ht)o<t<T, satisfying Assumptions (AC) and (EM). The process H has pre- 
dictable characteristics (0, C, v T *) or local characteristics (0, c, F *), and its 
canonical decomposition is 

H = J VcTdWj* + I J x(n H - v T *)(ds,dx), (3.1) 

OR 

where W * is a Pt„ -standard Brownian motion, fj, is the random measure 
associated with the jumps of H and v T * is the Py, -compensator of fi H . We 
further assume that the following conditions are in force. 
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(LR1): For any maturity T, there exists a bounded, continuous, deter- 
ministic function A(-, Ti) : [0, Ti] — > M, which represents the volatility 
of the forward LIBOR rate process L(-, T). Moreover, 

N 

5^|A(a,r0|<M, 
i=i 

for all s G [0,T*], where M is the constant from Assumption (EM), 
and X(s, T) = for all s > T,. 
(LR2): The initial term structure B(0,T), 1 < i < iV + 1, is strictly 
positive and strictly decreasing. Consequently, the initial term struc- 
ture of forward LIBOR rates is given, for 1 < i < N, by 

The construction of the Levy LIBOR model starts by postulating that 
the dynamics of the forward LIBOR rate with the longest maturity L(-,Tn) 
is driven by the time-inhomogeneous Levy process H and evolve as a mar- 
tingale under the terminal forward measure Py, . Then, the dynamics of the 
LIBOR rates for the preceding maturities are constructed by backward in- 
duction; they are driven by the same process H and evolve as martingales 
under their corresponding forward measures. 

Let us denote by Py <+1 the forward measure associated to the settlement 
date Tj + i, i G {0, . . . , iV}. The dynamics of the forward LIBOR rate L(-, Tj), 
for an arbitrary T, is given by 

L(t, Ti) = L(0, Ti) exp I f b L (s, Ti)ds + f X(s, T t )dHf +1 J , (3.3) 

where H i+1 is a special semimartingale with canonical decomposition 
t t 

H? i+1 = f ^F s dwT 1+1 + f f x{n H - v T ^)(ds, dx). (3.4) 

R 

Here W Ti+1 is a Pt-,, -standard Brownian motion and u Ti+1 is the Ft.,- 
compensator of \i H . The dynamics of an arbitrary LIBOR rate again evolves 
as a martingale under its corresponding forward measure; therefore, we spec- 
ify the drift term of the forward LIBOR process L(-,T) as 

b L (s,T i ) = -l\ 2 (8,Tjc. - J (e x ^ x - 1 - X(s,Ti)x)FT i+1 (dx). (3.5) 

R 

The forward measure Pr, +1 , which is defined on (Cl, T, (J r t)o<t<T i+1 ), is 
related to the terminal forward measure Pj 1 , via 

(3.6) 
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The Pr i+1 -Brownian motion W Ti+1 is related to the P^-Brownian motion 



via 



/ 



Wf + 1 = Wj w - j a(s, T i+l )^/F s ds = ... 
o 



where 



The PT i+1 -compensator of // , 1/ <+1 , is related to the Py, -compensator of 
/x via 

z/ Ti+1 (ds,dx) = /3(s,x,T i+ i)i/ T!+2 (ds,dx) = ... 

= ( f[ /3( S ,x,T0]^(d S ,dx), (3.9) 



\l=i+l 



where 



Remark 3.1. Notice that the process H i+1 , driving the forward LIBOR 
rate L(-,Ti), and H = H * have the same martingale part and differ only 
in the finite variation part (drift). An application of Girsanov's theorem for 
semimartingales yields that the P^ ^ -finite variation part of H is 



1 v T *(ds,dx). 



fc s J2 a{s,T l )ds+ j fxi f[ P{8,x,Ti)-l\ 

l=i+1 K V/=i+l / 

3.2. Option pricing and computational problems. The main scope 
of a market model for interest rate derivatives is to adequately describe 
the dynamics of interest rates as they are reflected in prices of derivatives. 
Hence, a good market model should be easily calibrated to option prices of 
liquid derivatives, i.e. caps and at-the-money swaptions. Calibration requires 
the fast computation of option prices; either in closed- form or using semi- 
analytical methods (e.g. Fourier transform methods). 

However, herein lies a major pitfall of the Levy LIBOR model: the process 
H Ti+1 driving the dynamics of L(-, Tj) is not a Levy process under Pr !+1 , or 
any other forward measure. One just has to observe that the compensator 
u Ti + 1 is random and not deterministic. Therefore, the characteristic function 
of the random variable H t I+1 is not available and Fourier methods cannot be 
used for option pricing. In other words, on top of the well-known problems 
of LMMs in pricing swaptions and other multi-LIBOR products, when the 
driving process has jumps even caplets cannot be priced in closed or semi- 
analytic form. 
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A remedy has been proposed in Eberlein and Ozkan (2005) and further 
refined by Kluge (2005). They propose to "freeze" the random terms in 
the compensator, i.e. to replace them by their deterministic initial values. 
The approximate process is then a Levy process and Fourier methods for 
option pricing can be applied. This method is equivalent to the "frozen drift" 
approximation, which does not however yield acceptable results. 

3.3. Terminal measure dynamics. Once closed-form or semi-analytical 
methods are not available for option pricing, a Monte Carlo simulation is 
the next alternative. In this section we derive the dynamics of LIBOR rates 
under the terminal measure. This is the appropriate measure for simulations 
in LMMs. 

Starting with the dynamics of the LIBOR rate L(-,Ti) under the forward 
martingale measure Pr i+1 , and using the connection between the forward 
and terminal martingale measures (cf. eqs. (3.7)-(3.10) and Remark 3.1), 
we have that the dynamics of the LIBOR rate L(-,Tj) under the terminal 
measure is given by 

L(t, Ti) = L(0, Ti) exp f b(s, T,)ds + f \{s, T t )dH s , (3.11) 

where H = (i?t)o<t<T, is the Pj^-time-inhomogeneous Levy process driving 
the LIBOR rates, cf. (3.1). The drift term 6(-,Tj) has the form 

b(s,T l ) = -]-\ 2 (s,T l )c s -c s X(s,T l ) J2 /' L }n ,Tl L X ^ T 

. f((eW> x -l) f[ /3(s,x,T0-A(s,T^xW*(dx), 

(3.12) 

where f3(s,x,Ti) is given by (3.10). Note that the drift term of (3.11) is 
random, therefore the log-LIBOR is a general semimartingale, and not a 
Levy process. Of course, L(-,Ti) is not a Pt„ -martingale, unless i = N 
(where we use the conventions Yl,i=i = an d lL=i = !)■ 

The equation for the dynamics under the terminal measure contains the 
next numerical problem, well-known for LMMs. The drift 6(-,Tj) depends 
on all subsequent LIBOR rates in the tenor, yielding a dependence structure 
that has the form of a triangular matrix; cf. Table 3.1. In other words, in 
order to simulate any rate one has to start by simulating the last rate, save 
the path and proceed iteratively. This means that simulations are very slow, 
while the burden on the random access memory (RAM) is also significant. 

The standard remedy to this problem is the so-called '"frozen drift" ap- 
proximation, where one replaces the random terms in the drift by their de- 
terministic initial values. This simplifies the simulations considerably, since 
rates are no longer state-dependent and can be simulated in parallel. How- 
ever, this approximation is very crude and does not yield acceptable results; 
cf. section 6. 
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L(t,?\) L(t,T 2 ) ... L(t,T k ) ... L(t,T N - 2 ) L{t,T N --,) L{t,T N ) 

L{t,T N ) L{t,T N ) ... L(t,T N ) ... L(t,T N ) L{t,T N ) 

L(t,T N -i) L(t,T N -i) ... L(t,7V_i) ... L(t,Tjv_i) 

L(t,T N - 2 ) L(t,T N - 2 ) 



L(t,T k+1 ) 

L(t,T a ) L(t,T s ) ... 

L(t,T 2 ) 

Table 3.1. Matrix of dependencies for LIBOR rates 



Moreover, an additional numerical problem arises in LMMs with jumps 
from the product term Y\i /?(-, -, T{). This term grows exponentially as a func- 
tion of tenor dates N and makes the simulations even more time-consuming. 

4. PlCARD APPROXIMATION AND DRIFT EXPANSION 

for the Levy LIBOR model 

The aim of this section is to derive approximation schemes for LIBOR 
models that can overcome the pitfalls of the the model - namely, the slow 
Monte Carlo simulations and the exponential growth of the product term. 
Firstly, using the idea of Picard iterations for the solution of SDEs, we 
derive approximate equations for the dynamics of LIBOR rates which are 
suitable for parallel computing. Secondly, by expanding and truncating the 
product term we can reduce the exponential to quadratic growth. Numerical 
examples show that these methods yield significant gain in computational 
time, while the loss in accuracy is very small. 

4.1. Picard iterations. In order to derive approximation schemes for the 
dynamics of LIBOR rates it is more convenient to work with the logarithm 
of rates. Let us denote by Z the log-LIBOR rates, that is 

Z(t,Ti) := log L(t,Ti) 

t t 



Z(0, Ti) + J b(s, Ti)ds + i A(s, T t )dH s , (4.1) 





where Z(0,Tj) = logL(0,T;) for all i G {l,...,N}. We can immediately 
deduce that Z(-,Tj) is a semimartingale and its triplet of predictable char- 
acteristics under P T ,, T(Z(-,Ti)|P T J = (5\C\i/), is described by 



B i = f b(s,Ti)ds 

Jo 



'0 

&= [ \ 2 ( S ,Ti)c s ds (4.2) 

Jo 

l A (x) * z/ = l A (A(s, Ti)x) *v T % A G B(R\{0}). 

The assertion follows from the canonical decomposition of a semimartingale 
and the triplet of characteristics of the stochastic integral process; see, for 
example, Proposition 1.3 in Papapantoleon (2007). 
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Remark 4.1. Note that the martingale part of Z(-,T{), i.e. the stochastic 
integral J Q X(s, Ti)dH s , is a time-inhomogeneous Levy process. However, the 
random drift term destroys the Levy property of Z(-,Tj), as the increments 
are no longer independent. 

The dynamics of log-LIBOR rates can be alternatively described as the 
solution to the following linear SDE 

dZ(t, Ti) = b(t, T z ; Z(t))dt + X(t, Ti)dH t 
Z(0,T l )=logL(0,T i ) (43) 

for all i E {1, . . . , iV} and all t £ [0, Tj\. We have introduced the term Z{-) in 
the drift term b{-, 2$; Z(-)) to make explicit that the log-LIBOR rates depend 
on all subsequent rates in the tenor. 

The idea behind the Picard approximation scheme is to approximate the 
dynamics of LIBOR rates by the Picard iterations for the SDE (4.3). The 
first Picard iteration for (4.3) is simply the initial value, i.e. 

Z^(t,T l ) = Z(0,T i ), (4.4) 

while the second Picard iteration is 

t t 

zW(t,T i ) = Z(0,T l )+ fb(s,T i ;Z^(s))ds+ I X(s,Ti)dH s 



I t 

= Z(0, Ti) + I b(s, Ti- Z(0))ds + I X(s, Ti)dH s . (4.5) 



Since the drift term b{-,Ti\ ^(0)) is deterministic, as the random terms have 
been replaced with their initial values, we can easily deduce that the second 
Picard iterate Z^'(-,Ti) is a Levy process. 

Remark 4.2. Comparing (4.5) with (4.1) it becomes evident that we are ap- 
proximating the semimartingale Z(-,Tj) with the time-inhomogeneous Levy 
process Z^ l '(-,Ti). 

4.2. Application to LIBOR models. We will now use the Picard iter- 
ations in order to deduce strong - i.e. pathwise - approximation schemes 
for the dynamics of LIBOR rates. More specifically, we will use the Picard 
iterates as proxies for the log-LIBOR rates in the drift term of the dynamics, 
cf. (3.12). Obviously, using the first Picard iterate Z' ' in the drift term we 
have just recovered the "frozen drift" approximation. 

Let us denote by Z{-,Ti) the approximate log-LIBOR rate stemming from 
using the second Picard iterate Z^ l > . The dynamics of the approximate log- 
LIBOR rate is 

t t 

Z(t,Ti) = Z(0,T<) + jb{s 1 T i -Z^{s))ds+ j X(s,Ti)dH s , (4.6) 
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where the drift term is provided by 

b( S ,T l ;Z^(s)) = -l\ 2 (s,T t )c s -c s X(s,T l ) jh - *'f ' f \{s,T x ) 

N \ 

V(-^>-l) n £(*,*, TO-AfoTOs Uf (dx), 

(4.7) 
with 

£(*,*, T„) = ^^W-'^)) ( e A(t,T,)x _ A + L (48) 

1 } l + 5 l exp(ZW(t-,T l ))\ J V ; 

The main advantage of this Picard approximation is that the resulting 
SDE for Z(-,Ti) can be simulated more easily than the equation for Z(-,Tj). 
Indeed, contrary to the dynamics of Z(-,Ti), the dynamics of Z(-, Tj) depend 
only on the Levy processes Z^'(-,Ti), i+1 < I < N, which are independent of 
each other. Compare also the "dependence matrix" for the approximate rates 
(Table 4.1) with Table 3.1. Hence, we can use parallel computing to simulate 
all approximate LIBOR rates simultaneously. This significantly increases 
the speed of the Monte Carlo simulations while, as the numerical examples 
reveal, the empirical performance is very satisfactory. 

Lft,Ti) L{t,T 2 ) ... L(t,T k ) ... L^.Tjy-i) L(t,T N ) 

Z^(t,T N ) Z^(t,T N ) ... Z^(t,T N ) ... Z^(t,T N ) 

Z^(t,T N ^) zW(t,Ttf-i) ... Z<V{t,T N -i) ... 

zU(t,T N - 2 ) zW(t,T N - 2 ) 

Z (1 \t,T k+1 ) 

zV>(t,T„) Z^(t,T 3 ) ... 

Z w (t,T 2 ) 
Table 4.1. Matrix of dependencies for approximate LIBOR rates 



Remark 4.3. Note that the Picard approximation can be also used in 
case one wants to apply P(I)DE methods for the valuation of derivatives 
in LIBOR models, and yields an analogous simplification of the problem. 

Remark 4.4. Let us point out that the Picard approximation scheme (4.6)- 
(4.8) we have developed is universal and can be applied to any LMM. We 
can replace the Levy process H driving the dynamics of LIBOR rates by 
a general semimartingale X with random predictable characteristics (thus 
also incorporating stochastic volatility). Subject to certain assumptions X 
has the canonical decomposition 

t t 

X t = J ^F s dW s + i I ' x(n H - i/)(da,dx); (4.9) 

OR 
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compare with (3.4). Then we can construct an LMM driven by the semi- 
martingale X following the steps for the Levy LMM in Eberlein and Ozkan 
(2005). We can also follow analogously all the steps for the Picard approxi- 
mation; in this case, the Picard iterate Z^- 1 ' will have the same dynamics as 
in (4.5), with H replaced by X. 

4.3. Drift expansion. This part is devoted to the integral term in the drift 
of the LIBOR dynamics; see (3.12) again. Obviously this is a problem solely 
related to LMMs driven by jump processes. 

Let us introduce the following shorthand notation for convenience: 

A,:=A(s,T,) and U := L(s,T t ). (4.10) 

We denote by A the part of the drift term that stems from the jumps, i.e. 

A = /((e--l) fl (^^(e^-^-A.xW). 

R V l=i+1 ' 

(4.11) 

In theory, one could simply employ a straightforward numerical integration 
to compute A. However this is not feasible in practice since a numerical 
integration should be performed at each step of the Monte Carlo simulation. 
An alternative solution is to express A in terms of the cumulant generating 
function of (the jump part) of H. 

Observe that a product of the form J^J i=1 (l + ol{) appears, where ol\ := 
1+5 L ( e ^ lX ~ !)• This product can be expressed in terms of so-called elemen- 
tary symmetric polynomials. Let k < N, then the elementary symmetric 
polynomial of degree k in N variables is given by 

e fe (ai,..., aw) = Yl a h x ■■■ x a ik . (4.12) 

l<h<-<i k <N 

Hence we have that 

N 

~[(1 + ai) = l + £i(ai,...,cw) H \- e N (ai, . . . ,a N ). (4.13) 

1=1 

One can immediately deduce that, while the drift term stemming from 
the diffusion part is a first order polynomial in 1 jg l L , A is an iV-th order 
polynomial. More importantly, the number of terms on the RHS of (4.13) 
is 2 N . Hence, we need to perform 2^ computations in order to calculate 
the drift of the LIBOR rates. Since N is the length of the tenor, it becomes 
apparent that this calculation is feasible only for short tenors. If, for example, 
N = 40 this amounts to more than 1 trillion computations. 

In order to make this computation more feasible we will truncate the RHS 
of (4.13) at the first or second order. The first order approximation of the 
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product term is 

A!= [ ((e XiX - l) (l + ei(a m , ...,a N )- X iX )) F s T *(dx) 



v^ s i L i („\,x 



Ds, 

= / (e A ' x - 1 - Xix) Fj" (dx) 

+ E T ffi;/( ei ' I - 1 )(^- 1 )^^ 

1=1+1 R 

TV . ,. 

= K W+ E 7^tt( k ( a * + a ')- k M-^)), (414) 

i=i+l ' ' 

and the order of the error is 

A = A' + 0(iV||L|| 2 ). (4.15) 

Similarly the second order approximation is provided by 

A" = J ((e A '* - l) (l + (ei+e 2 )(a l+1 ,. ..,a N )- A,x)) if*(dx) 

R 

N 

' + SiLi 



;(Ai)+ Yl ! + ^ L , ( K ( A * + A + K ( A + K ( A * 

Z=i+1 



51 r 



_ + SiLi 1 + (5j.-Lt 



+ 

i+l<fe 

x (k(A, + A/ + Afc) - n(Xi + Aj) - k(A» + Afc) 

- K (A fe + A/) + «(Ai) + k(Aj) + «(A fc )) , (4.16) 

and the order of the error is 

A = A" + 0(iV 2 ||L|| 3 ). (4.17) 

Since the LIBOR rate is an order of magnitude smaller than the number 
of tenor dates, these approximations are justified. Indeed, numerical results 
show that truncation at the second order yields very satisfying results, while 
the gain in computational time is very significant. 

5. Derivatives on LIBOR rates 

In this section we briefly describe the derivatives we will use for the numer- 
ical illustration. Namely we use caplets and swaptions, which are the most 
liquid derivatives in the interest rate markets. We will also use forward rate 
agreements (FRAs) as a benchmark for the different approximations because 
they have model independent values. Of course, since we have developed a 
pathwise approximation, we could also consider many other derivatives - 
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especially path-dependent ones - for the illustration. We avoid this for the 
sake of brevity. 

5.1. Caplets. The price of a caplet with strike K maturing at time T, using 
the relationship between the terminal and the forward measures (3.6), can 
be expressed as 

Cq(K, Ti) = 5iB(0, T i+l ) E PTi+i \{L(T u Ti) - K)+] 

rdPT 
= 5B(0,T l+1 )TE FTf [^^\^(L(T u T l )-Ky 

N 

= <5B(0,T»)]Ep T .[ n (l + SL(Ti,Ti))(HTi,T^-K)+^. (5.1) 

This equation will provide the actual prices of caplets corresponding to sim- 
ulating the full SDE (Euler discretization) for the LIBOR rates. In order to 
calculate the Picard approximation prices for a caplet we have to replace 
L(-,T.) in (5.1) with L(-,T). 

5.2. Swaptions. A payer (resp. receiver) swaption can be viewed as a put 
(resp. call) option on a coupon bond with exercise price 1; cf. section 16.2.3 
and 16.3.2 in Musiela and Rutkowski (1997). Consider a payer swaption 
with strike rate K, where the underlying swap starts at time T and matures 
at T m (i < m < N). The time-Tj value is 

(m \ + 

1- ^ c k B{Ti,T k )\ 
k=i+l J 

/ m / k— 1 

= '>- E(«n r 



i 



(5.2) 



Ck 



(5.3) 



-+sm, Tl) 

where 

SfrK, i + l<k<m — 1, 

1 + SfrK, k = m. 

Then, the time-0 value of the swaption is obtained by taking the P^ 
expectation of its time-Tj value, that is 

S = S (K,Ti,T m ) 

' m / k—1 



B(0,Ti)Ep 
B(0, %) 

x E Pt>> 



k=i+l 



1 



;=•; 



SL(Ti,T t ) 



N f m / fc-1 

n(i+«w.n» (i- E {^UTTmf, 

l=i \ k=i+l v l=i v l 



Ti) 



hence 



S = 5(0,T,)IEp Tt 



m y N s ' 



(5.4) 
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where c% := —1. 

5.3. Forward Rate Agreements. A forward rate agreement with strike 
K and notional value of 1 with expiry at time Tj, has the value of 5{{K — 
L(Ti,Ti)) at expiry. The time zero value of the contract has the model in- 
dependent value of 

¥ (K, Ti) = 5 l B(0, T i+l )[K - L(0, Tj)], (5.5) 

or zero if the contract is struck at-the-money (K = L(0,T)). In the fol- 
lowing section we will compare the known true values with simulated prices 
generated from the terminal measure expectation of the payoff, i.e. 

N 

¥ (K, Ti) = SB(0, T*) E Ptj [ J] (l + SL(T t , T,)) (K - L(T U Tj))] . (5.6) 

6. Numerical illustration 

The aim of this section is to demonstrate the accuracy and efficiency of 
the Picard approximation scheme and the drift expansions in the valuation 
of options in the Levy LIBOR model. In the first section we demonstrate 
the accuracy of our methods compared to a standard Euler discretization 
of the LIBOR SDE in pricing caplets and swaptions. In the second section 
we estimate the speed of our method; finally we compare with alternative 
approximations. 

We will consider a simple example with a flat volatility structure of 
A(-,T) = 18% and zero coupon rates generated from a flat term structure 
of interest rates: 5(0, Tj) = exp(— 0.04 • Tj). We consider a tenor structure 
with 6 month increments (Si = g). 

The driving Levy process H is a normal inverse Gaussian (NIG) process 
with parameters a = 5 = 12 and (j, = /3 = 0, resulting in a process with 
mean zero and variance 1. We denote by \i the random measure of jumps 
of H and by u(dt,dx) = F(dx)dt the Pt, -compensator of fi , where F is 
the Levy measure of the NIG process. The necessary conditions are then 
satisfied for term structures up to 30 years of length because M = a, hence 
E£i|A(-,rO| = 10.8<a. 

The NIG Levy process is a pure-jump process with canonical decomposi- 
tion 

H = j I x(fi H - i/)(da, dx). (6.1) 

R 

The cumulant generating function of the NIG distribution, for all u E C 
with |5ftu| < a, is 



k(u) = 5a — 5\/a 2 — u 2 . (6-2) 

6.1. Accuracy of the methods. The Picard approximation should be con- 
sidered primarily as a parallelizable alternative to the standard Euler dis- 
cretization of the model. This will therefore be the benchmark to which we 
compare. In order to avoid Monte Carlo error we use the same discretization 
grid (5 steps per tenor increment) and the same pseudo random numbers 
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(50000 paths) for each method. The pseudo random numbers are gener- 
ated from the NIG distribution using the standard methodology described 
in Glasserman (2003). 

Starting with caplets we can see in Figure 6.1 the difference between the 
Euler discretization and the Picard approximation expressed in price (left) 
and implied volatility (right). The difference in price is small with max. er- 
rors around half of a percentage of a basis point. On the right we see some- 
what larger errors with a maximum slightly below 1 basis point of implied 
volatility for low strike mid-maturity caplets. Implied volatility is normally 
quoted in units of 1 basis point while bid-ask spreads are usually around 
at least 5 bp of implied volatility. The errors are therefore at acceptable 
levels. Note also that in experiments not shown we found that the levels and 
patterns of the errors are insensitive to the number of discretization points 
as well the number of paths. 

The errors display a non-monotonic behavior as a function of maturity 
with peaks around mid-maturity. The non-monotonicity can be explained 
by the fact that volatility dominates the price of options in the short end, 
making the drift, and any error in it, less relevant. As maturity increases 
the importance of the drift grows relative to volatility but the state depen- 
dence becomes less critical as the number of "live" rates decreases. These 
two opposing effects result in the mid-maturity peak that we observe. This 
pattern is also noted in the study by Joshi and Stacey (2008) . 

As we established in (4.13) the number of terms needed to calculate the 
drift grows at a rate 2 N . In market applications N is often as high as 60 
reflecting a 30 year term structure with a 6 month tenor increment. At 
this level even the calculation of one drift term becomes infeasible and this 
necessitates the use of the approximations introduced in (4.14) and (4.16). 
We investigate the errors introduced by the drift expansion by comparing 
them with the full numerical solution obtained as before using the true drift 
from (3.12). The results are plotted in Figure 6.2. Here we can see that 
the first order drift expansion adds errors of fairly low magnitude, whereas 
the second order drift expansion performs significantly better with errors so 
small that they can be disregarded. We also plot the Picard approximation 
alone as well as combined with the second order drift expansion and these 
are similarly indistinguishable. 

Continuing on with swaptions, in Figure 6.3 we plot again differences in 
implied volatility between our methods and the Euler discretization. We 
observe even smaller errors than for caplets, most likely because swaptions 
span a broad range of maturities which smooths out the higher errors in the 
mid-maturity region. 

6.2. Computational speed. In terms of computational time the largest 
gain by far is realized when using the drift expansions in (4.14) and (4.16). 
In the example above the CPU time for the full numerical solution is 1.5 
hours; after applying the first order and second order drift expansion it drops 
to 1.3 and 27.2 seconds respectively. 

The Picard approximation in itself does not contribute to the compu- 
tational speed unless parallelization is employed. In fact, it is slightly but 
insignificantly slower as the auxiliary processes Z^ 1 ' have to be evolved along 
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Figure 6.1. Difference in price (left) and implied caplet 
volatility (right) between the Euler discretization and the 
Picard approximation (in basis points). 
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Figure 6.2. Difference in implied caplet volatility (in basis 
points) between Euler discretization and 4 other methods. 



with the rates. On the left in Figure 6.4, we plot CPU time as a function 
of the number of paths for the Picard approximation and the full Euler dis- 
cretization. In both cases the second order drift approximation scheme in 
(4.16) is employed. The computations are done in Matlab running on an 
Intel i7 processor with the capability of running 8 processes simultaneously. 
Here we see the typical linear behavior as the number of paths are increased; 
notice though that the Picard approximation has significantly lower slope. 
On the right in Figure 6.4 we plot CPU time as a function of the number of 
rates. One can see CPU time exponentially increasing, revealing that large 
gains in computational time are realizable when using the Picard approxi- 
mation scheme and the drift expansion. 

Needless to say, the speed advantages of the Picard approximation are 
only partially revealed in these graphs since we are using desktop computer; 
further speed increases can of course be realized as the access to more CPUs 
(or clusters of PCs) becomes available. This is already part of the infras- 
tructure of many large financial institutions. 
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Figure 6.3. Difference in implied swaption volatility (basis 
points) between Euler discretization and 4 other methods. 
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Figure 6.4. CPU time as a function of the number of paths 
(left) and the number of rates N (right). 



6.3. Comparison with other methods. Unfortunately, very little work 
has been done in the area of approximations for LMMs driven by general 
semimartingales, leaving us without any standard method to compare with, 
other than the frozen drift approximation. As mentioned in the introduction 
the existing work has focused mainly on the log-normal case and the case 
of finite intensity jump-diffusion models. However, some of the techniques 
applied to the log-normal case can be adapted to our setup as well. 
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Assume we want to simulate the LIBOR rates from time t to time t-\- h, 
where t + h < Ti . We have 

/ t+h t+h \ 

L(t + h, Ti) = L(t, Ti) exp f b(s, T t ; Z(s))ds + f X(s, T t )dH s , 



I 



where 6(-,Tj; Z(-)) is the state dependent drift function defined in (3.12); cf. 
also (4.3). The standard (log)-Euler scheme leads to 

t+h 

b(s, Tf, Z(s))ds « b(t, T t ; Z(t)) h. (6.3) 

t 
This can be further refined as noted first in Kloeden and Platen (1999) (see 
also Hunter, Jackel, and Joshi 2001, and Joshi and Stacey 2008 and the ref- 
erences therein) by using instead 

t+h 

I b(s, T t ; Z(s))ds « (h(t, T % - Z(t)) + h(t + h, Ti- Z(t + h))\ h. (6.4) 

t 
The second term in the parenthesis requires the knowledge of the LIBOR 
rates L(t + /i,Tj + i), . . . , L(t + h, T/v); one therefore has to simulate these 
rates. This can be done in a separate simulation step and the procedure 
is known in the LMM literature as predictor- corrector (PC) method. One 
can also note that when rate i is evolved under the terminal measure it only 
depends on rates k > i. Furthermore, if we start with i = N we have no state- 
dependence in the drift. We can then generate realizations of L(t + h, Tjv) 
without discretization error and these can be used in the drift of rate iV — 1 
as described above. Realizations of rate N — 1 can then be generated with 
the corrected drift from (3.12), which can be subsequently used in the drift 
of rate N — 2, and so forth. This latter method is referred to as iterative 
predictor- corrector (IPC). 

IPC has been found to often outperform PC in the log-normal case studied 
in Joshi and Stacey (2008). It is also slightly more efficient since it does not 
require a separate simulation step for the rates at time t + h. 

Remark 6.1. One should point out that PC and IPC are alternative dis- 
cretization schemes to the Euler discretization. 

We can also combine PC and IPC with the Picard approximation by 
merely using b(-,Tf, Z^\-)) instead of b(-,Tf,Z(-)) in (6.4). Furthermore, 
PC and IPC will actually be equivalent when applied to the Picard ap- 
proximation since the drift term does not involve the rates themselves, but 
the auxiliary processes Z^- 1 ', which makes the order in which the rates are 
evolved irrelevant (see section 4.2). 

We compare PC and IPC with our methods in 3 different cases. Since 
we do not have the true price of caplets or swaptions we instead compare 
prices of at-the-money forward rate agreements (FRA) since these all have 
a known model independent price of zero. Here we keep Monte Carlo error 
sufficiently low by simulating 5 million paths, employing antithetic sampling 
as the only variance redaction measure. Looking at the top left of Figure 
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6.5 we have kept the discretization grid dense at 5 steps per tenor period. 
Here we see prices which for all methods are sufficiently small at max levels 
of 2.5% of a basis point. Furthermore the different methods are more or less 
indistinguishable - certainly in statistical terms with 95% confidence limit 
halfwidths ranging from about 0.004 bp in the short end to to 0.03 basis 
points in the long end. The frozen drift case is left out in these graphs since 
the errors/prices are so big that they dominate all other methods. 

The second graph in the top right of Figure 6.5 shows prices for a dis- 
cretization grid of 1 step per tenor period. Here the prices for Picard and 
Euler clearly reflect the higher discretization error and the PC and IPC 
methods are indistinguishable, and significantly lower than Picard and Eu- 
ler. Note that prices are around the same level as in the 5 step per tenor 
period discretization. 

Finally, we price each FRA using a single long discretization step from 
time zero to the maturity of each contract, also referred to as "long-stepping" . 
In this case the Picard and Euler methods are analogous to the frozen drift 
approximation, while PC and Picard+PC are equivalent to each other; hence 
we draw only PC, IPC, and frozen drift. We see, somewhat surprisingly that 
the frozen drift slightly outperforms PC and IPC which contradicts the re- 
sults in the log-normal case previously studied in the literature. Nevertheless 
the errors are quite high and beyond an acceptable level for all methods. 

The general conclusion one can draw from these graphs is that Picard+PC 
with 1 discretization step per tenor period would be the preferable choice. 
The errors are indistinguishable from regular PC and IPC, but the method 
has the advantage that prices can be parallelized in the maturity dimension, 
and the gains in computational time shown in the previous section can be 
realized. 



7. Conclusion 

This paper derives new approximation methods for Monte Carlo simula- 
tion in LIBOR models. The methods address the problem of speed in Monte 
Carlo simulations by allowing for parallel computing through Picard ap- 
proximations. In particular, our method decouples the interdependence of 
the rates when moving them forward in time in a simulation, meaning that 
computations can be parallelized in the maturity dimension. Furthermore, 
the largest computational load in the model stems from the exponential 
growth of the drift terms. We reduce this growth to quadratic through trun- 
cated expansions of the product terms that appear in the drift. The accuracy 
is very high if the second order expansion is employed and we showed that 
it reduces the computational load immensely. 

Numerical methods for LIBOR market models driven by general semi- 
martingales are still in their infancy. As we have demonstrated there is still 
work to be done in this area, in particular developing algorithms for pricing 
derivatives using long time-stepping. Predictor-corrector methods do not 
perform very well in this case and a finer discretization grid needs to be 
employed. 
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